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PREFACE 


The three-dimensional rigid-lid model is intended to be used for 
hydrothermal predictions of closed basins subjected to a heated discharge 
together with various other Inflows and outflows. This volume has been 
written in order to assist any prospective user In applying the model to 
specific sites. Derivation of the governing equations and various other 
details have been omitted. The programs are fairly general and only 
one subroutine and a data file has to be rewritten for specific cases. 

This work was sponsored by the National Aeronautics and Space 
Administration (NASA-KSC) and the Environmental Protection Aqency 
(EPA-RTP) . 
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ABSTRACT 


The three-dimensional rigid-lid model was developed by the thermal 
pollution group at the University of Miami and verified for accuracy at 
various sites. The model results have been found to be fairly accurate 
in all the verification runs. The model is intended to be used as a 
predictive tool in future sites and this manual has been written to enable 
any user to be able to apply it without difficulty. 
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SECTION 1 


INTRODUCTION 


The need for mathematical modeling in predicting and monitoring 
thermal pollution was discussed in previous reports by Veziroglu et al. 
(1973, 1974). Predictive studies of ecosystems can only be made by 
mathematical models. A prior knowledge of the effects of disturbances 
is essential for environmental impact studies. Thus, the mathematical 
model is a crucial tool In decisions Involving power plant siting, land 
development, etc. 

The University of Miami team undertook development of a methodology 
using remote sensing and numerical modeling to study thermal pollution. 
The use of remotely- sensed data in modeling has been discussed by Sen- 
gupta et al. (1974). The remote sensing effort has been discussed in . 
detail in previous publications. This volume has been written so as to 
enable a user to apply the mathematical model to new sites for predictive 
purposes. 

The hydrodynamics and thermodynamics of an ecosystem are con- 
trolled by geometry, meteorological conditions and physical characteristics 
of the water such as density, salinity and turbidity. In this model the 
effects of salinity and turbidity have been neglected. Hence, the govern- 
ing equations are composed of the three-dimensional Navier-Stokes equa- 
tions and the energy equation. Various assumptions can be made for 
different situations leading to simplification or elimination of equations. 

The main simplifying assumption In this case is the rigid-lid assumption. 
This means that surface height fluctuations are not simulated by this 
model, and this is a reasonable assumption for most applications (e.g.. 
Lakes) . 

The rigid-lid model has the following capabilities: 

1. It predicts the wind-driven circulation. 

2. It predicts the circulation caused by inflows and outflows to the 
domain. 

3. It predicts the thermal effects in the domain. 

4. It combines the aforementioned processes. 

The calibration procedure consists of comparing ground-truth cor- 
rected airborne radiometer data with surface isotherms predicted by the 
model. 
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SECTION 2 


RECOMMENDATIONS 


Various numerical models have been developed to study the effects 
of heated discharge and meteorological conditions on bodies of water. 

Most of these models are one or two dimensional. These models have a 
high computational speed but only give horizontally or vertically averaged 
values of temperatures. 

Three-dimensional models, however, have a much finer resolution 
but they consume larger computer time. The three-dimensional rigid- 
lid model can be used to obtain detailed temperature and velocity distri- 
butions in a domain where surface gravity waves are small compared to 
the depth of the domain. This model, as compared to free-surface mo- 
dels, runs faster since surface gravity waves are eliminated by the 
rigid-lid assumption. 

A proper method of using this model would be to run a one-dimen- 
sional model initially to obtain a rough picture of the temperatures and 
then using this model to obtain a better resolution, the 1-D results being 
used as ambient conditions. 

The following improvements have been suggested for the model. 

1. Since all natural flows are turbulent, proper turbulent closures are 
needed to make the model meaningful. At present, the simplest 
possible closures, namely constant eddy viscosities and eddy diffu- 
sivities, have been used. However, better results may be obtained 
by using a higher order closure. 

2. At present, the model uses uniform horizontal grids and stretched 
vertical grids. Nonuniform horizontal grids could be introduced for 
better resolution near the boundaries. 

3. The program has been written to be run as a batch-job on the com- 
puter. It could be made interactive so as to enable the user to run 
it on a terminal. However, this would require some modifications in 
order to reduce the storage space. 
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SECTION 3 


PROCRAM DESCRIPTION AND FLOW CHART 


DESCRIPTION OF PROCRAM ALCORITHM 

The governing equations for a body of water which are derived from 
the basic laws of conservation of mass, momentum and energy are shown 
in Table 1. These equations incorporate a vertical ly- stretched coordinate 
system so as to make the model general enough to handle any kind of 
bottom topography. The problem is set up as an initial value problem. 
The initial values of the water velocities and temperatures are specified 
and the model is run so as to give the values of the above quantities in 
subsequent time periods using an explicit scheme. The sequence of the 
calculations are as follows: 

1. The initial values of the velocities and temperatures are read into 
the program, the region of interest within the basin being classified 
into interior, corner or boundary points. (Subroutines used are 
READ 3K, INITIA, INITIT, HEIGHT.) 

2. The data, which includes the boundary conditions such as the 
various meteorological parameters like surface wind speed, air 
temperature, humidity and solar radiation are read into the program 
using subroutine READ 2. 

3. Depending on the site chosen, the various discharges (volume flow 
rate, velocities and temperatures) in and out of the basin are read 
into the model. These are incorporated in the subroutine INLET1. 

4. The momentum, continuity and energy equations are now solved to 

determine the velocities and temperatures in the subsequent time 
steps. The predictive equation for pressure (viz., the Poisson 
equation) is solved iteratively to determine the pressures at various 
points of the domain. (Note: Because of the rigid-lid assumption, 

the surface or lid pressure is no longer atmospheric.) 

THE PROGRAM FLOW CHART IS SHOWN IN FIGURE 1 

The various subroutines used are as well as a brief description of 
their functions are shown in Tables 2 and 3. 


Symbols Used in Governing Equations 
(Quantities with bar are dimensional) 


p = density 
f = temperature 

1 = t(u !t + v fr’ * 
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* /U ref 



V = 

* /U ref 
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Table 1. 


Governing Equations 


Continuity Equation: 


3(hu) 3(hv) . 3J3 ... 

3a 3 B 3 y 


Momentum Equation: 


3(hu) 3(huu) . 3(huv) , u 3(fJu) h 

~TT ~ + "TF n 3 y " R b 


3 P 

= -hrr - hB : + r 


1 3/1 3 A x 1 3 flJV. 

3^ h “J + R 7af h TT ) 

e e 


and 



3(hv) 3(huv) 3(hvv) h 3 ( 3v) h 

3t 3 a 33 3y R B 





* 



if h l 


3 v, 


E*R^ h 3 y V3 y ’ 
e 

Hydrostatic Equation: 

fy * E u (,+ » )h 

Energy Equation: 

3 (hT) 3 (huT) 3 (hvT) . 3 (AT) 

3 t 3 a 3 3 ^Ty 


- 1 2_/h 3T i 

* ~ 7af h TS’ 

e 


3_/ h 3T, 

F# h 7i ? 


P~& 

e 


1 

R 


3 y v3 Y 
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SECTION 4 


LIST OF PROGRAM SYMBOLS USED IN MAIN PROGRAM 
DESCRIPTION OF MAIN VARIABLES 

A. A - constant In equation of state, p = A + BT + CT J 
AREF - reference eddy viscosity 

AA - value of 'V 1 at plume inlet 

ABR - 1 /Rossby number 

AH - 1 /Reynolds number 

Al - coefficient In front of pressure term 

AKT - (K s ) CH ref ) /CB z ) 

AP - coefficient in front of pressure term 

ARBP - arbitrary pressure 

AV - — t 4 — where e = r 

c^Rg L 

A3 - normalized vertical eddy coefficient of viscosity 
ANGLE - wind direction angle 

B. B - constant in equation of state, o - A + BT + CT 2 
BB * value of 'V' at plume inlet (at 1=10) 

B2 - °C p B v 

BV - normalized vertical eddy diffusivity, normalized with respect 
to reference eddy diffusivity 

C. C - constant in equation of state, p = A + BT + CT 2 

CC - value of y (constant) 

CW - temperature gradient at vertical boundaries 

CB - temperature gradient at the bottom 
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D - U 

at previous 

time step 

D1TZ - 

3T 

TT 


DPX - 

3 P 
3 x 

• 

DRY - 

3 P 
3y 


DPSX - 

3X 


DPSY - 

3P s 

ay 


DT - time increment 

DX - increment in 

^’direction 

DY - increment in 

y-direction 

DZ - increment in 

Z-direction 

D1HUX 

3 (hu) 
3 x 


D1HVY 

- 3 (hv) 



3y 


D1HUUX 

3(huu) 
3 x 


D1HUVY 

3(huv) 

3y 


D1HVVY 

3(hvv) 

3y 


D1UY - 

3 U 

3y 


D1VX - 

3 v 
3 x 


D2UX - 

3 2 u 
3 X 3 


D2VX - 

3 2 v 
3 x 1 


D1VWX 

3(vw) 

' TT 


D1UZ - 

3 u 

TT 
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D2UZ 


3*u 

TzT 


D1VZ 
D2VZ 
D1A3Z 
DLZ - 


3 v 
3 z 

3 J v 
TTT 




E. E - 
EPS 
EUL 
EX - 

F. FH - 
FW - 

C. C - 

H. H - 
HI - 
HREF 
HX - 

HY - 

I. IN - 
IWN ■ 


V at previous time step 

- convergence criterion 

- Enler number 

residual error in pressure iteration 

forcing function in pressure equation 

factor in wind stress calculation formula 

dummy variable for V (for future time step) 

dummy variable for U (for future time step) 

h 

nondimensional depth = p 

- reference depth 

3 H 

3 a 

3 H 
3 6 

maximum number of grid points in x-dlrectlon 

- maximum number of half-grid points in x-direction, IWN = 
IN - 1 


I - Index of x-axis, main grid 
ITN - Index for number of Iterations 


IW - index for x-axis, half grid 

IRUN - index for number of runs 
* 0, first run 

= 1, from second time onwards 
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ISCNX, ISGNY - determine signs of TAUX and TAUX respectively 

J. J - index for y-axis, main grid 
JW - index for y-axis, half grid 

JWN - maximum number of half-grid points in y-direction 
JWN - JN - 1 

JN - maximum number of main grid points in y-direction 

K. K - Index for Z-axis 

KSTORE - specified usage of tape for storing results 

KN - maximum number of main grid points in Z-direction 

KISS - surface heat transfer coefficient (nondimensional) 

L - maximum length of the domain 

LN - number of time steps to be computed 

LLN - total number of time steps /LN 

M. MAR - number to describe general location of a point in the main 
grid 

MRH - number of describe general location of a point in the half 
grid 

MAXIT - maximum number of iterations 

O. OMEGA - relaxation factor 

P. P - nondimensional pressure 

PN - New pressure, nondimensional 

PI NTH - dummy variable for pressure (future time step) 

R. R - dimensional density at main grid points 
RE - Reynolds number 

RB - Rossby number 

RINTX - density integrated with respect to x 
RINTY - density integrated with respect to y 
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RO - nondimensional density at main grid points 
ROW - nondimensional density at half grid points 
RREF - reference density (gm/cc) 

RW - dimensional density at half grid points (gm/cc) 

RADN - solar radiation (w/m J ) 

T. T - nondimensional temperature at main grid points 
TO - initial temperature (dimensional) (°C) 

TAMB - ambient temperature (dimensional) (°C) 

TAIR - air temperature (dimensional) (°C) 

TAI - coefficient in front of convective terms in the energy equa 
tion, = 1. 

TAH - J- where P = R x 

r _ e e r 

e 

TAV - p"^~T where e = £ 

e 

TE - equilibrium temperature (dimensional) (°C) 

TTOT - total time elapsed 
TAUX - 3 u / 3 y (nondimensional) 

TAUY - 3v/3y (nondimensional) 

TEM - cfimensional temperature at main grid points 
TEMW - dimensional temperature at half-grid points 
TREF - reference temperature 

TW - nondimensional temperature at half-grid points 
ILL - temperature at the discharge point (nondimensional) 

TSU - water surface temperature (nondimensional) 

TDEW - dewpoint temperature (dimensional) 

U. U - velocity in x-dlrection (nondimensional) 

V. V - velocity In y-direction (nondimensional) 
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VVIS - vertical eddy viscosity (nondimensional) 

W. W - velocity in Z-dlrection (nondimensional) 

WH - W at half-grid points 

WHLDT - time derivative of WH at lid (i.e., y^(WH) II - 0) 

X. XINT - integral of x terms on the right-hand side of Poisson's 

equation 

X - horizontal coordinate across discharge 

Y. YINT - integral of y terms on the right-hand side of Poisson's 

equation 

Y - horizontal coordinate across discharge 

Z. Z - vertical coordinate 
MARKER MATRICES 

The following number convention is used for the MAR = matrix sys- 
tem, which classifies points (or nodes) on the main grid system = 

(Refer to Figure). 

MAR = 0. points outside the region of interest. 

MAR = 1, point on the far y-boundary. 

MAR = 2, point on the near y-boundary. 

MAR = 3, point on the near x-boundary. 

MAR - 4, point on the far x-boundary. 

MAR = 5, outside corner on near x-boundary and far y-boundary. 

MAR = 6, inside corner on far x-boundary and far y-boundary. 

MAR = 7, outside corner on near x-boundary and near y-boundary. 

MAR = 8, inside corner on near x-boundary and near y-boundary. 

MAR = 9, outside corner on far x-boundary and near y-boundary. 

MAR = 10, outside corner on far x-boundary and far y-boundary. 

MAR = 11, points in the interior of the region of interest. 
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trix 

MRH 

MRH 

MRH 

MRH 

MRH 

MRH 

MRH 


The following number convention is used to describe the MRH (ma- 
for the half-grid system) . 

= 1, corner at far x-bovindary and far y-boundary. 

- 2, points on near y-boundary. 

= 3, points on near x-boundary. 

= 4, corner at near x and near y-boundaries. 

= 6, far corner on x-axis. 

= 7, corner at far x and y-boundarles. 

= 9, interior grid points. 
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SECTION S 


PREPARATION OF RUNS 


This section presents the steps to be followed in order to run the 

model for a particular location. 

1. The boundaries are chosen depending on the particular situation, the 
general idea being to include all inflows and outflows. If a heated 
discharge enters the body of water the region of interest must be 
chosen so as to include this since it is a major factor in determining 
the size and spread of the resulting plume. 

2. The grid size is chosen depending on the resolution required. The 
user should remember that the choice of the grid size directly deter- 
mines the maximum allowable time step since this is directly re- 
lated by the various stability criteria. (See choice of time step in 
Section 6.) 

3. Specify number of full-grid points IN, JN, KN and number of half- 
grid points IWN, JWN. Since the actual domain may be smaller than 
the total rectangular region, INxJNxKN, the marker matrices MAR 
and MRH are used to specify the domain so that points outside the 
domain of interest skip the subsequent calculations. 

4. IRUN is specified (- 0 for the first run, = 1 for subsequent runs). 
KSTORE is specified to indicate whether any tape has been assigned 
to store results of the run. 

KSTORE = 0 if no tape has been assigned. 

= 1 if tape has been assigned. 

LLN is specified to denote the number of hours of simulation to be 
carried out. 

5. The depths at various places within the domain are specified using 
subroutine HEIGHT. The various Inflows and outflows to the domain 
are specified using INLET 1. (For details please refer to Biscayne 
Bay run, Sengupta et al. (1975).) 

6. The various data like solar radiation, wind speed, wind direction and 
dewpoint temperature are specified in a data file which is made by 
the main program. 

For further details see the next section. 


13 



SECTION 6 


INPUT DATA 


The data that Is required for the execution of the main program Is 
listed In Table 3 In the order It appears. Note, the data Input symbols 
have already been defined in Section 4. Moreover, the following remarks 
should be observed. 

* Free format Is used for all data input. 

* Distinction must be made for integer and real number. 

* The order of the cards must be followed. 
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SECTION 7 


PLOTTING PROGRAMS 


The plotting programs for the 3-D rlgld-lld model are distinct from 
the main program and subroutines used to run it. The user has an 
option of either using a tape (Unit 8) during running the main program 
TMAINN to store the results or just run it without storing the results. 

For making subsequent continuation runs of TMAINN all that is required 
is the result of the last hour in the previous run. For plotting, however, 
one needs the results of ail the hours for which results are to be plotted. 
These results are used as input data to run the various plotting programs. 

DESCRIPTION OF PLOT PROGRAMS 


The following are the main plotting programs. 
PLOT - plots surface isotherms. 

PLUV - plots u, v components of the velocities (i.e. , 

PLUW - plots u, w components of the velocities (i.e., 

PLVW - plots v, w components of the velocities (I.e., 

SUBROUTINES 


'K' sections). 
'I* sections). 
*i* sections). 


The various plot programs and subroutines are shown in Table 4. 

Other subroutines seen in these programs (e.g., ARROHD, FLINE, 
etc.) are standard FORTRAN subroutines used for plotting, using a 
CALCOMP x,y plotter, and are hence ommltted in the above listing. 
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APPENDIX A 


EXAMPLE CASE 


INTRODUCTION 

The area of Interest is Lake Keowee in South Carolina, which was 
formed from 1 968 through 1971 by damming the Little and Keowee rivers. 
The lake is located about 40 km west of Creenvitle and constitutes Duke 
Power Company's Keowee-Toxaway complex. 

Lake Keowee has two arms connected by a canal (maximum depth 
30.5 m) . There are three power plants on the lake, namely, the Oconee 
Nuclear Station, Keowee hydro station and Jocassee-pumped storage 
station. The Oconee Nuclear Station is a three unit steam-electric sta- 
tion with an installed capacity of generating 2580 MW. The Oconee 
Nuclear Station draws in condenser-cooling water from the lower arm 
of Lake Keowee and discharges the heated effluent to the upper arm of 
the lake. The intake structure for the condenser-cooling water allows 
water from 20 to 27 m depth (full pond) to pass through. The discharge 
structure has an opening from 9 to 12 meters below the water surface 
(full porid) through which the CCW returns directly to the upper branch 
of the lake. 

Lake Jocassee is located north of Lake Keowee and is used as a reser- 
voir for Jocassee-pumped storage station. Lake Keowee also serves as 
the lower pond for this station. The Jocassee station has reversible 
turbines with a maximum generating flow (into Lake Keowee) of about 
820 m 3 / sec and a maximum pumping flow (out of Lake Keowee into Lake 
Jocassee) of about 775 m 3 /sec, the net flow into Lake Keowee from Jocas- 
see being about 15.5 m 3 /sec. 

Lake Keowee has a full pond elevation of 243. 8 m above MSL. At 
full pond it has a volume of approximately 1.18 x 10 y m 3 , an area of 
74 km 2 , a mean depth of 15. 8 m and a shoreline of about 480 km. The 
outflow from Lake Keowee is through Keowee hydro station and may vary 
from approximately 1.4 m 3 /sec (leakage) to 560 m 3 /sec. Maximum allow- 
able draw-down of the lake Is 7.6 m. 

A map of the area of Interest Is shown in Figure 3. 

PROBLEM STATEMENT 

The objective of the present work is to find the three-dimensional 
temperature and velocity distributions in the region where the effects 
of the thermal discharge are noticeable. The effects of Jocassee-pumped 
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storage station, Keowee hydro station as well as the meteorological condi- 
tions have been incorporated. 

The region of interest is chosen to include the effects of the Oconej 
Nuclear Station discharge, the outflow through Keowee dam and the impact 
of the J ocas see- pumped storage station on the velocity and temperature 
distributions in Lake Keowee. The depth of the domain is cut off at 1G 
meters, since this is the level at which the thermocline occurs. Hence, 
for running the model, a constant depth region is considered. The plan 
view of the domain is shown in Figure 4. (Note: For variable depth 

refer to Biscayne Bay simulation studies by the University of Miami 
thermal pollution group.) In this figure, AB is an open boundary which 
takes care of the flow from or to the Jocassee-pumped storage station. 

'C' shows the position of the flow in the canal connecting the two arms 
of the lake. 'D' is the discharge point for the Oconee Nuclear Station 
and 'E' is the outflow from Keowee hydro station. 

The inclusion of the above results in a domain 2895.6 m x 2438.4 m 
in the horizontal plane. The horizontal grid size (in x and y directions) 

Is 152. 4 m x 152.4 m. giving a total of 20 x 17 (= 340) nodes in the 
horizontal plane, out of which 293 lie in the region of interest. The 
16 m constant depth region of interest is divided into 4 equal slices of 

4 m each, giving at total of 5 nodes in the vertical (Z) /direction. 

Hence, there are 293 x 5 nodes (grid points) in the region of interest. 
This region is specified using the MAR and MRH marker matrices (Figure 

5 and Figure 6) . 

Boundary Conditions 

On the Jocassee effect boundary, the flow, .velocity (varying with 
time) is specified. Open-boundary condition (|— ■ = 0) is specified for 
the temperature. y 

The same is done for the Keowee hydro boundary. The only differ- 
ence is that the values specified are at three points in the vertical plane 
(i.e., at K = 1, 2 and 3) since this region covers the discharge area. 

For the Oconee Nuclear Station, the discharge velocity as well as 
the discharge temperature is specified at the discharge point. 

Open-boundary conditions are specified for the temperature and 
velocity at the canal. This, however, leads to a possible violation of 
mass balance in the region of interest. This mass unbalance will 
actually show up as a variation in the water level in the lake which is 
beyond the capability of the rigid-lid model. 

At all solid boundaries as well as the artificial bottom (since the 
bottom Is cut off at 16 m) perfect insulation (temperature gradient = 0) 
and zero velocity conditions are assumed. 

At the surface, the vertical component of the velocity is specified 
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as zero (rigid-lid constraint) . Surface wind shear stress and heat 
transfer coefficient are specified. 

Initial Conditions 

The initial values of the water velocities are assumed to be zero. 
The initial temperature of the lake is assumed to be equal to the am- 
bient water temperature (determined by running a one-dimensional 
model) and is taken to be uniform throughout the domain. 

CALCULATION OF PARAMETERS AND INPUT DATA 

Reference Quantities 

Reference length = L = maximum length of the domain = 2895.6 m. 
Reference horizontal eddy viscosity A . = 0.002 L 4/3 

= 38311.48 cm 2 sec. 

For better agreement with data the value chosen is 60,000 cm 2 /sec. 
Reference depth = H = 16 m. 

Reference vertical A v = 0.002 x (H) 4/3. 

Eddy viscosity = 37. 43 cm 2 /sec. 

Reference velocity = V ref = 30 cm /sec. 

Reference temperature = T re ^ = 10.0°C. 

Reference time = L/V re ^ = 9652 sec. 

Calculation of Inflows and Outflows into the Domain (Used in INLET!) 


Oconee Nuclear Station Discharge Velocity — 

The discharge is considered to take place through a point at a 
depth of 12 m (k = 3). The discharge velocity is calculated as follows: 


“ 

12 m 

> 

f 




152. 4 m 


The total discharge into the basin is equal to: 

(100 2-2 x V x 152.4 x 12) = Q 
m 
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where Q = average discharge In m 3 /sec 


.'.V 


8144. 1 
60 


m / sec 


= 7. 42207 cm /sec 

The average values of Q over 24 hrs Is taken since the variation 
is negligible. 

V V 

Nondimensional discharge velocity = y — = -jj = 0. 24740 


Keowee Hydro Discharge Velocity — 

The outflow through the Keowee hydro station is through a channel 
1 52. 4 m x 12 m. 

The volume flowrate Q = (152. 4 x 12 x V) m 3 /sec 
where V = discharge velocity (m/sec) 

. * . V = (Q / (1 52. 4 x 12)] m/sec = ( i S2 . 2x T 6o ) cm/sec 

Q is specified as a function of time In INLET1. 

The procedure for nondimensional! zation is similar. 

Jocassee Flow Velocity — 

The entire flow to or from the Jocassee-pumped storage station is 
assumed to take place through the entire upper boundary (AB in Figure 
4). The flow through this area (shown below) Is assumed to be uniform 
and is assumed to take place simultaneously with the outflow through the 
Jocassee station. 


16 m 


V = Q /[ (16x1 3x152. 4) X 100] cm/sec. 

Q = flow through Jocassee (m 3 /sec). 

Q is positive when Jocassee Is generating (i.e., the flow is into 
the region of interest) and negative when pumping (i.e., flow out of 
region of interest) . 



13 x 152. 4 m 
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SAMPLE INPUT 


The following are the Inputs to TMAINN contained in the data file 
IPUT (which includes values calculated earlier). 


Input 


# 

1 


2 


3 


No. of Data 
In Card 


Symbol 


3 IRUN 
KSTORE 

LLN 

2 VVIS 

ABR 

4 Al 
AH 

AV 


Value 


= 0 
= 1 
= 3 


= 37.43/60,000 = 0.00062384 

= 0. 78 

= 1.0 


60,000 


= ( 


30 x 2895 x 100 
2895. 6 2 


= 0.01228172 


16 


) AH = 402. 08304 


AP 

4 4 EPS 

MAXIT 

OMECA 

ARBP 

5 3 DX 

DY 

DZ 

6 3 TAI 

TAH 

TAV 

7 3 A 

B 


= 1.0 
- 0.001 
= 60 
= 1.8 
= 1.0 

= 1 52. 4/2895.6 = 0.05263 

= 0.05263 

= 4/16 = 0. 25 

= 1.0 

= AH = 0.01228172 
= AV = 402.08304 
= 1.000428 

= -0.000019 


t 
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Input No. of Data Cw „. . 
# In Card Symbo1 

c 


8 

1 

TO 

9 

3 

EUL 



CW 



CB 

10 

2 

AA 



CC 

n 

1 

TLL 

12 

1 

TAU 

13 

1 

DT 


1 4 1 CTTOT 

15 1 ISOTOP 

16 6 WS 

TSU 

TDEW 

RADN 

ISCNX 

ISGNY 

17 1 ANCLE 
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LAKE KEOWEE APPLICATION-EXECUTION DECK 

The following execution deck Is for use in the UNIVAC 1100 computer 
at the University of Miami. These may have to be modified if a different 
computer is used. 

(ALL PROGRAMS AND SUBPROCRAMS COMPILED AND STORED IN FILE) 
First Run 

1. 9 ASG, AX FILE. 

(THE FILE IS ASSICNED FOR THE RUNI 

2. @ ASG,T 8, 1 6N , TAPENAME. 

(A TAPE FILE NAMES '8' IS BEING ASSIGNED. THE TAPE IS 
9-TRACK, AND THE REEL NUMBER IS 'TAPENAME') 

3. 0 PRT,S FILE. TMAINN 

(THE MAIN PROGRAM IS PRINTED) 

4. 0 PACK FILE. 

(THE FILE IS PACKED) 

5. 0 PREP FILE. 

(ENTRY POINT TABLE IS PREPARED) 

6. 0 MAP.S 

7. IN FILE. TMAINN 

8. LIB FILE. 

9. END 

10. 0 XQT 

11. 0 

(VALUE FOR IRUN, FIRST RUN: IRUN=0) 

12. 24 

(NUMBER OF HOURS REQUIRED, MINIMUM=1 HOUR, MAX=24) 
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13 . 0 


(0 IF MAGNETIC TAPE IS REQUIRED TO STORE RESULT. IF 
NOT, ANY NUMBER) 


14. 0 ADD FILE. INPUT 

(INPUT DATA FILE FOR THE PARTICULAR RUN) 

15. 0 FIN 

EXECUTION DECK FOR PLOT PROGRAMS 

1. 0 ASC.AX FILE. 

2. 0 ASC.T 8., 16N, TAPENAME. 

3. 0 ASC.T 11., 16N, PLOTTAPE. 

(A MACNETIC TAPE FILE NAMED '11' IS BEING ASSIGNED. 
TAPE IS 7-TRACK AND THE REEL NUMBER IS 'PLOTTAPE'. 
PLOTS ARE STORED ON THIS TAPE) 

4. 0 PRT.S FILE. PLOTTER 

(THE PLOT PROGRAM IS PRINTED) 

5. 0 PACK FILE. 

6. 0 PREP FILE. 

7. 0MAP.S 

8. IN FILE. PLOTTER 

9. LIB FILE. 

10. END 

11. 0 XQT 

12. 0 ADD FILE. INPUT 

13. 0 FIN 


THE 

THE 
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Table 2. Subroutines Required in Main Program TMAINN 


No. 

Name 

Description 

Remarks 

1 

DVISV 

Computes D1VY, D2VY, D1VX 
and D2VX . 

Called by subroutine INTE. 
Schemes used similar to 
DVISU. 

2 

DVISU 

Computes D1UX, D2UX, and 
D1UY. 

Called by INTE. p, p 
are computed at interior, 
boundary or corner pts 
by scheme similar to the 
one used in DINERU. 

3 

DVVY 

Computes D1HVVY. 

Called by INTE. y- (hvv) 
is computed for interior, 
boundary or corner by a 
scheme similar to the one 
used in DINERU. 

4 

DUVY 

Computes D1HUVY. 

Called by INTE. y- (huv) 
is computed for interior, 
boundary and corner pts 
by a scheme similar to the 
one used in DINERU. 

5 

DINERU 

Computes D1HUUX and 
D1HUVX. 

Called by INTE. The re- 
sults are used in Poisson 
equation for pressure. 

6 

TPRINK 

Prints temperatures at a grid 
point. 

Called by TMAINN. 

7 

PRUV 

Prints the values of U and V 
at all main grid points. 

Called by TMAINN. 

8 

PRITEX 

Prints the No. of iterations 
(ITN) and final residual error 
in solving the Poisson equation 

Called by TMAINN. 

9 

TPRIN1 

Prints the input parameters. 

Called by TMAINN. 

10 

STORE 2 

Stores values of Input para- 
meters and physical quantities 
on tape #8 

Called by TMAINN. 

11 

RWR 

Computes real vertical veloci- 
ties from modified vertical 
velocities used in equations at 
integral grid points. 

Called by TMAINN. 
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Table 2. Subroutines Required in Main Program TMAINN (Continued) 


No. Name 


Description 


Remarks 


Computes real vertical veloci- 
ties at half-grid points. 


Called by TMAINN. 


DENSTY I Uses the equation of state and Called by TMAINN. 
computes density field from the 
temperature field. 

Allows for vertical mixing at a If the temp at the grid 
particular grid point. Program pt just above it is less 
is called by TMAINN. and the difference Is more 

than a specified maximum, 
the two temperatures are 
averaged. 

Sets the values of temperature 
field at time step 'n' equal to 
the temperature field at (n+1) 
after all computations for time 
step 'n 1 are completed. 


TEMB2 


TEMU 


Computes temperatures at the 
boundary points in the domain 
of interest. 


Called by TMAINN. 


Computes temperatures at the Called by TMAINN. 

interior points of the domain of 

interest. 


Computes vertical velocities at Called by TMAINN. 
half-grid points. 


OLDUV Sets the values of D and E 

equal to U and V respectively 
in order to retain values of U 
and V at one time step lag. 

UVTOP Computes U and V at the top 
using wind stress boundary 
conditions. 


Called by TMAINN. 


Called by TMAINN. Com- 
putations are made for 
MAR = 11 only (internal 
grid points). 


Computes U and V for variable Called by TMAINN. 
density at successive time 
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Table 2. Subroutines Required in Main Program TMAINN (Continued) 


No. 

Name 

Description 

Remarks 

22 

PRE1L 

Computes pressure for far 
field from Poisson’s Equation 
at half-grid points. 

Called by TMAINN. 

23 

FORCE 

Computes R.H.S. of Poisson's 
Equation at half-grid points. 

Called by TMAINN. 

24 

DPSXY 

Computes DPSX and DPSY. 

Called by TMAINN. 

25 

ROINTY 

Computes Y in the Poisson's 
Equation. p 

Called by TMAINN. 

26 

ROINTX 

Computes X in the Poisson's 
Equation. p 

Called by TMAINN. 

27 

CORINT 

Adds integral of Coriolis' 
component XI NT and YINT. 

Called by TMAINN. 

28 

INTE 

Computes XINT, YINT, DPSX, 
and DPSY. 

Called by INTE. 

29 

WHATIJ 

Computes the values of W at 
1, J from the values of WH at 
IW, JW. 

Called by TMAINN. 

30 

WHTOP 

Sets the value of WH equal to 
zero at the surface. 

Called by TMAINN. 

31 

ERROR 

Calculates "Hlrt and Harlow" 
correction term at half-grid 
points and at the surface 
(WHLDT) . 

Called by TMAINN. 

32 

READ2 

Reads in input parameters and 
physical quantities stored on 
tape #7. 

Corresponds to store 2. 
Called in by TMAINN. 

33 

INLET 1 

Puts in velocities u and v 
pheme discharge, etc. into 
the model. 

Called by TMAINN. 

34 

HEIGHT 

Inputs depths of the basin 
into the model. 

This subroutine is for a 
constant depths model. 
Called by TMAINN. 
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Table 2. Subroutines Required In Main Program TMAINN (Continued) 


No. 

Name 

Description 

Remarks 

35 

INITIT 

Sets initial temperature field. 

Sets the temperature field 
equal to ref temp at all 
grid points. 

Called by TMAINN. 

36 

INITIA 

Initializes values of U, V, WH, 
W, D, E and PI NTH. 

Called by TMAINN. 

37 

READ 3K 

Classifies region of interest 
into interior, corner and 
boundary points using matrix 
MAR. 

Called by TMAINN. 

38 

IPUT 

Data files containing values of 
input data for the respective 
days. 
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Read data cards 


First run? (IRUN = 0} 


READ 2 - Reads tape for 
data of previous run 


I N LET 1 - Inputs V, U, 
TD to start run 


READ3K - Classifies the re- 
gion into interior, corner or 
boundary points. 

I N IT I A - Initializes U, V, W, 
WH, D, E and P. 

I N IT IT - Initializes T and p. 
HEICHT - Inputs depth. 
INLET1 - Inputs V, U, TD 
TD to start run. 


ERROR 

_ 

Computes WHLDT 

WHTOP 

- 

Sets WH equal to zero 

WHATIJ 

- 

Computes W 

INTE 

- 

Computes XINT, YINT, DPSX and DPSY 

ROINTX 

- 

Computes X p 

ROINTY 

- 

Computes Yn 

DPSXY 

- 

Computes DPSX, DPSY 

FORCE 

- 

Computes FH 

PRE1L 

- 

Computes PN 

UVTOP 

- 

Computes U and V at the top 

OLDUV 

- 

Sets D = U and E = V 

RWH 

- 

Computes vertical velocities at half-grid points 

WHATIJ 

- 

Computes W at grid points 

TEMI4 


Computes temperatures at interior points 

TEMB2 

- 

Computes temperatures at boundary points 

OLDT 

- 

Updates temperatures 

TEQB 

- 

Allows for vertical mixing 

DENSTY 

- 

Computes density field 

INLET1 


Inputs V, U and TD every time step 


A- new- hour? 


, Yes 


RWRH - Computes real vertical velocities 

RWR - Computes real vertical velocities 

STORE2 - Stores computed results on magnetic tape (unit 8) 

TPRIN1 - Prints input parameters 

PRITEX - Prints number of iterations and error 


Figure 1. Flow chart (main program) 
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Figure 1 (Continued) . Flow chart (main program) 
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Figure 2. Coordinate and grid system 
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Table 3. Input Data to TMAINN 


Input 

# 

No. of Data 
In Card 

Symbol 

. ■ — - 

Definition /Value 

1 

3 

1RUN 

= 

0 for first run 



LLN 

s 

No of hours of simulation 



KSTORE 

= 

0 if no tape is assigned 

1 if tape is assigned 

2 

2 

VVIS 


Nondimensional vertical eddy 
viscosity 

3 

4 

ABR 

Al 

= 

1 /Rossby No. = n— — 
ref 

Coefficient in front of Inertia 
term =1.0 



AH 

— 

1 /Reynolds No. = 

Ref eddy hoz viscosity 

tt-ttl 

ref 


4 

AV 

= 

(1/e J Re) (c = H/L) 



AP 

= 

Coefficient in fron of pressure 
term =1.0 

4 

4 

EPS 

= 

Convergence factor = 0.001 



MAXIT 

2 

Maximum number of Iterations 
for Poisson Equation 



OMEGA 

s 

Relaxation factor = 1.8 



ARBP 

= 

Arbitrary pressure = 1.0 

5 

3 

DX 

s 

Horizontal grid spacing (x dir.) 



DY 

2 

Horizontal grid spacing (y dir.) 
Ay/L 



DZ 

s 

Vertical grid spacing (z dir.) 
Az/H 

6 

3 

TAI 

2 

Coefficient of convective terms 
in energy equation = 1.0 



TAH 

= 

Horizontal eddy diffusivlty 
AH (usually) 
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Table 3. Input Data to TMAINN (Continued) 


Input 

# 

No. of Data 
In Card 

Symbol 

Definition/Value 



TAV 

- 

Vertical eddy diffuslvlty 




55 

AV (usually) 

7 

3 

A 

- 

1.000428 These are coefficients 




S 

-0.000019 in the equation of 




= 

-0.0000046 state for water where 
p = A + BT + CT 2 (gm/cc) 

8 

1 

TO 

= 

Reference temperature (°C) 

9 

3 

EUL 

= 

Euler No. = v j 

1 rer 



CW 

= 

Temperature gradient at vertical 
boundary 



CB 

= 

Temperature gradient at the 
bottom 

10 

2 

AA 


Nondimensional discharge velocity 



= 

(discharge velocity) /U ref 



CC 

- 

No dimensional depth = h/H ref 

11 

1 

TLL 


Nondimensional discharge tem- 
perature = (T d - T q ) fX Q 

12 

1 

TAU 


Surface shear stress (from Wilson 
Curve) (Refer to Figure 7) 

13 

1 

DT 

= 

Nondimensional time step 



= 

iT<L/U re( ) 

14 

1 

CTTOT 

s 

Converts nondimensional time to 
hours 

15 

1 

ISTOP 

s 

Number of hours of previous run 

16 

6 

WS 

= 

Wind speed (m/sec) 



TSU 

s 

Air temperature (°C) 



TDEW 

= 

Dewpoint temperature (°C) 



RADN 

s 

Incident solar radiation (w/m J ) 
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Table 3. Input Data to TMAINN (Continued) 


Input 

# 

No. of Data 
In Card 

Symbol 

Definition /Value 



ISGNX 

= +1 if x component of W is 
negative s 

= -1 If x component of W Is 
positive s 

= +1 If y component of W s is 
negative 

= -1 If y component of is 

positive 

17 

1 

ANGLE 

= Direction of W (degrees) with 
respect to the s x axis 
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Table 4. Plotting Programs 


No. 

Name 

Program Description 

Remarks 

1 

PLOT 

Plots surface Isotherms 


2 

PLUV 

Plots velocities, K section 


3 

PLUW 

Plots velocities, j section 


4 

PLVW 

Plots velocities, i section 


5 

ECHKON 

Calculates equal temperature points 

Called by PLOT 

6 

CONLIN 

Draws the isotherms 

Called by ECHKON 

7 

ENDER 

Writes the values of the tempera- 
ture on the isotherms 

Called by ECHKON 
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Figure 3. Lake Keowee 





Figure *». Lake Keowee (region of interest) showing inputs and outputs (for 3-D model) 



Figure 5. MAR marker matrix 
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Figure 6. MRH marker matrix 
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Table 5. Meteorological Data for Lake Keowee (February 27, 1979) 


Time 

(hrs from 
midnight) 

Wind Speed 
(m/s) 

Air Temp 
(°C) 

Dewpoint 

Temp 

(°C) 

Solar 

Radiation 

(w/m 2 ) 

Wind 

Direction 

(Degrees) 

1 

1.833 

-0. 33 

-2.78 

0.0 

15° 

2 

1.073 

-0.72 

-1.67 

0.0 

75° 

3 

2. 325 

-1.61 

-1.61 

0.0 

O 

o 

to 

4 

1. 565 

-2. 22 

-2. 28 

0.0 

15° 

5 

2.056 

-1.83 

-1.89 

0.0 

50° 

6 

1.788 

-2.17 

-2. 22 

0.0 

85° 

7 

2.012 

-2. 72 

-2.78 

20.94 

85° 

8 

2. 280 

-1.67 

-2. 78 

1 95. 39 

60° 

9 

0.626 

0.01 

-3. 33 

369. 85 

5° 

10 

1. 386 

3.06 

-2. 22 

544. 31 

75° 

11 

1.609 

5. 83 

-2. 22 

655. 31 

15° 

12 

1.788 

8. 83 

-1.39 

725. 75 

40° 

13 

3.129 

11.06 

-2.78 

746.68 

80° 

14 

2. 593 

12. 28 

-5.00 

704.81 

70° 

15 

I. 520 

13. 39 

-5.56 

579. 20 

80° 

16 

1. 207 

1 3. 89 

-5. 56 

383. 81 

75° 

17 

1.565 

13. 83 

-5.61 

146.55 

55° 

18 

1. 609 

13.72 

-3. 33 

20. 94 

15° 

19 

2.056 

11.72 

-4. 44 

0.0 

30° 

20 

1.162 

9.72 

-2.78 

0.0 

25° 

21 

1.772 

8. 33 

5. 28 


55° 

22 

2. 861 

7.78 

5. 56 

0.0 

55° 

23 

2. 995 

7.00 

5. 28 


50° 

24 

1. 386 

5. 28 

3. 89 

0.0 

60° 
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Table 6. Inflows and Outflows to Lake Keowee 


Time 

Feb. 27, 1978 

Oconee 
Discharge 
(m J /min) 

Oconee 
Discharge 
Temp (°C) 

Net Jocassee 
Flow 

(C.F.S.) 

Keowee Hydro 
Flow 

(C.F.S.) 

12.00 a.m. 

7505. 3 

18. 6 

-14395 

48 

1.00 

7498. 1 

18. 5 

-18754 

48 

2.00 

7492.0 

18. 4 

-18805 

48 

3.00 

7492.0 

18. 5 

-18713 

48 

4.00 

7491.6 

18.3 

-18698 

48 

5.00 

7494. 3 

18. 3 

-18688 

48 

6.00 

7488. 2 

18. 3 

-15939 

48 

7.00 

7481. 8 

18.2 

3484 

3668 

8.00 

7485.6 

18. 3 

16823 

17540 

9.00 

7488. 2 

18. 2 

13503 

8488 

10.00 

7497.7 

18. 3 

5470 

8096 

11.00 

7504. 1 

18. 3 

| 100 

2680 

1 2. 00 p.m. 

7503. 4 

i8. 4 

100 

48 

1.00 

7506.0 

18. 5 

100 

48 

2.00 

7506. 4 

18. 5 

100 

48 

3.00 

7503. 4 

18. 4 

100 

48 

4.00 

7301. 9 

18. 4 

100 

48 

S.00 

7507. 5 

18. 4 

100 

48 

6.00 

7511.0 

18.4 

100 

48 

7.00 

7516. 2 

18.4 

100 

48 

8.00 

7518. 9 

18. 3 

100 

48 

9.00 

7520. 4 

18. 2 

100 

48 

10.00 

7516.6 

18. 2 

100 

48 

11.00 

7509. 4 

18. 2 

100 

48 

12.00 a.m. 

7507. 2 

18. 2 

-4382 

48 
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FLOW 


Figure 7. Keowee hydro discharge (February 27, 1979) 
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APPENDIX B 

FORTRAN SOURCE PROGRAM LISTING 
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LIST OF MAIN PROGRAM AND SUBROUTINES 
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SAMPLE PLOTS 
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ISOTHERMS RT K= 1. 

LAKE KEQWEE-( RIOIO-LIO MOOEU 
SIMULATIONS FOR FEB. 27 1979 
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ISOTHERMS AT K= 1. 

LAKE KE0WEE-( RIGIO-LIO MODEL) 
SIMULATIONS FOR FEB. 27 1979 


IDO 





RUN NO) L9Q 8. 

DISCHARGE VELOCITY * 7.42CM/3EC 
DISCHARGE TEMPERRTURE* 3l.7*C 
NINO SPEED (MAX) i 3.Q9M/SEC 
CURRENT! J0CAS3E FLOMJi 1.1 CM/SEC 
TOTAL SIMULATED TIME I 1.01 HRS 


LENGTH SCALE! METERS J 

0,-00 

01 .00 

VELOCITY SCALE! CM/SEC) 

0.00 

L 

12.00 




VELOCITIES RT K= 1 . 

LAKE KE0WEE-( RIGID-LID MODEL ) 
SIMULRTIQNS FOR FEB. 27 1979 


1 41 



RUN NO i L9Q 6. 

OISCHRRGE VEL8CITY i 7.42CM/SEC 
DISCHARGE TEMPERATURE i 31 .7®C 
NINO SPEED CMRX1 i 3.09M/SEC 
CURRENTC JQCRSSE FL0H)» 1.1 CM/SEC 
T9TRL SIMULATED TIME i 1.01 HRS 

* 10 * 

LENGTH SCRLEC METERS) ° .‘°° 6 . x ‘ 00 

VEL9CITY SCRLEC CM/SEC) ° .* Q0 L 2 * 00 




VELOCITIES AT K= 2. 

LAKE KEOWEE-t RIGID-LID MODEL ) 
SIMULATIONS FOR FEB. 27 1979 
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VELOCITIES AT K= 3. 

LAKE KE0MEE-( RIGID-LID MQOEL ) 
SIMULATIONS FOR FEB. 27 1979 


1 « 




RUN NQi L00 6. 

DISCHARGE VELSCITY t 7.42CM/SEC 
DISCHARGE TEMPERATURE i 31.7 C 
NINO SPEED (MAX) « 3.Q9M/SEC 
CURRENT( J0CA33E FL8H) » 1.1 CM/SEC 
TOTAL SIMULATED TIME * 1.01 HRS 



VELOCITIES AT K= 4. 

LAKE KE0WEE-I RIGID-LID MODEL ) 
SIMULATIONS F0R FEB. 27 1979 


TU4 



RUN NO i LQQ 8 . 

OISCHRROE VELOCITY * 7.42CM/SEC 
OISCHRROE TEIiPERRTURE i 31.7 C 
HIND SPEED CMRX) a 3.Q9M/SEC 
CURRENT! JQCASSE FL 8 HJ 1 1.1 CM/SEC 
T0TRL SIMULATED TIME a 1.01 HRS 



VELOCITIES AT K= 5. 

LAKE KEG WEE-( RIGID-LID M0DEL ) 
SIMULATIONS F0R FEB. 27 1979 
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